library(nlme)
library(lattice)
library(lindia)
library(mgcv)
library(tidyverse)
library(descr)
library(patchwork)
theme_set(theme_minimal())
pretty_estimates <- function(m) { # m: model object
cbind(Est. = m$coefficients, confint(m))
}
lang <- read_delim("http://publicifsv.sund.ku.dk/~lts/basal/data/varicer.txt", delim = " ") %>%
mutate(patient = factor(patient),
treat = factor(treat, levels = c("P", "A")))
lang5 <- lang %>%
filter(minutter <= 5)
ggplot(lang, aes(x = minutter, y = pressure, colour = patient)) +
geom_line(alpha = 0.5, show.legend = FALSE) +
geom_point(size = 1, show.legend = FALSE) +
facet_wrap(~ treat, nrow = 1)
lang %>%
group_by(treat) %>%
summarise(n_patients = n_distinct(patient))
## # A tibble: 2 x 2
## treat n_patients
## <fct> <int>
## 1 P 9
## 2 A 8
lang %>%
group_by(patient) %>%
summarise(n_obs = n())
## # A tibble: 17 x 2
## patient n_obs
## <fct> <int>
## 1 1 13
## 2 2 9
## 3 3 9
## 4 4 12
## 5 5 9
## 6 6 9
## 7 7 9
## 8 8 11
## 9 9 10
## 10 10 17
## 11 11 17
## 12 12 9
## 13 13 10
## 14 14 14
## 15 15 11
## 16 16 10
## 17 17 10
# Fewer and fewer patients, mainly in placebo arm (not a good sign)
lang %>%
group_by(minutter, treat) %>%
summarise(n_patients = n_distinct(patient)) %>%
ggplot(aes(x = minutter, y = n_patients, colour = treat)) +
geom_line() +
coord_cartesian(ylim = c(0, NA))
# Mean VP develops
lang %>%
group_by(treat, minutter) %>%
summarise(mean_vp = mean(pressure)) %>%
ggplot(aes(x = minutter, y = mean_vp, colour = treat)) +
geom_line() +
coord_cartesian(ylim = c(0, NA))
ggplot(lang, aes(x = minutter, y = pressure, colour = patient)) +
geom_line(alpha = 0.5, show.legend = FALSE) +
geom_point(size = 1, show.legend = FALSE) +
facet_wrap(~ treat, nrow = 1) +
# geom_label(aes(label = patient), ~ filter(., minutter == 0), size = 8/ggplot2::.pt) +
# facet_wrap(~ patient) +
coord_cartesian(xlim = c(0, 5)) # Q2
One patient is quite different from the rest in the active arm (suddent drop around 2 minutes in). To identify the person we can add labels to the lines (or plot them individually, see the two commented-out lines above) and we see it patient 13. Let’s remove this patient for now as it seems an outlier of some sort and difficult to model.
Restrict population and time span
lang <- lang %>%
filter(patient != 13, minutter <= 5)
# Plot mean over time by arm
plot_df <- lang %>%
group_by(treat, minutter) %>%
summarise(mean_vp = mean(pressure),
sd_vp = sd(pressure),
.groups = "drop")
plot_mean_vp <- ggplot(plot_df, aes(x = minutter, y = mean_vp, colour = treat)) +
geom_line() +
coord_cartesian(ylim = c(0, NA))
plot_sd_vp <- ggplot(plot_df, aes(x = minutter, y = sd_vp, colour = treat)) +
geom_line() +
coord_cartesian(ylim = c(0, NA))
plot_mean_vp /
plot_sd_vp
This is not expected: values drop too quickly after 3 minutes, due to the low patient count. Otherwise the standard deviation is pretty stable over timer and no too different between groups. This difference it likely that those in the placebo with higher presssures are lost.
At each follow-up time we lose patients, and because those lost patients have higher VP’s, the mean will of course decrease as well.
With respect to randomisation, we need to look at the baseline measurements, where the two groups seem quite comparable. Little differences are to be expected as it’s a little population and randomisation can’t account for random variance. We shouldn’t do a test per se, but if one were to do it it would be a T-test:
t.test(pressure ~ treat, data = filter(lang, minutter == 0))
##
## Welch Two Sample t-test
##
## data: pressure by treat
## t = -0.60929, df = 12.462, p-value = 0.5533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -9.687360 5.439741
## sample estimates:
## mean in group P mean in group A
## 19.23333 21.35714
We could test the difference at a specific time point, say, after 5 minutes:
t.test(pressure ~ treat, data = filter(lang, minutter == 5))
##
## Welch Two Sample t-test
##
## data: pressure by treat
## t = 0.10015, df = 4.4469, p-value = 0.9246
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -7.013492 7.560158
## sample estimates:
## mean in group P mean in group A
## 14.53333 14.26000
with(t.test(pressure ~ treat, data = filter(lang, minutter == 5)),
c(diff(rev(estimate)), conf.int))
## mean in group P
## 0.2733333 -7.0134915 7.5601582
table(filter(lang, minutter == 5)$treat) # very few observations, so non-parametric test perhaps better
##
## P A
## 3 5
wilcox.test(pressure ~ treat, data = filter(lang, minutter == 5))
##
## Wilcoxon rank sum test
##
## data: pressure by treat
## W = 9, p-value = 0.7857
## alternative hypothesis: true location shift is not equal to 0
Anyway, don’t read too much into this as we didn’t account to loss to follow-up. But we do start a table for the results:
result <- tibble(metode = "3d: T-test efter 5 minutter",
behandlingseffekt_5_minutter = "-0.27 (-7.57, 7.01)",
p_value = 0.92)
Two approaches:
lang <- lang %>%
mutate(dif05 = pressure - baseline)
t.test(dif05 ~ treat, data = filter(lang, minutter == 5))
##
## Welch Two Sample t-test
##
## data: dif05 by treat
## t = 3.8401, df = 4.4813, p-value = 0.01492
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.332268 12.881065
## sample estimates:
## mean in group P mean in group A
## 0.9666667 -6.6400000
result <- result %>%
rbind(c("3e: T-test på ændring 0-5 minutter", "-7.61 (-12.88, -2.33)", 0.015))
Almost nothing changes in placebo group, but rather large change in active arm. We add the signs to results to make them consistent.
ancova <- lm(pressure ~ baseline + treat, data = filter(lang, minutter == 5))
summary(ancova)
##
## Call:
## lm(formula = pressure ~ baseline + treat, data = filter(lang,
## minutter == 5))
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 0.2870 -0.0642 -0.2229 1.0531 4.1460 -5.4169 -1.6232 1.8411
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4249 3.6829 1.473 0.2007
## baseline 0.6714 0.2328 2.884 0.0344 *
## treatA -5.1968 2.9418 -1.767 0.1376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.28 on 5 degrees of freedom
## Multiple R-squared: 0.6249, Adjusted R-squared: 0.4748
## F-statistic: 4.164 on 2 and 5 DF, p-value: 0.08619
pretty_estimates(ancova)
## Est. 2.5 % 97.5 %
## (Intercept) 5.4249496 -4.04223058 14.892130
## baseline 0.6713796 0.07288985 1.269869
## treatA -5.1967840 -12.75881549 2.365248
result <- result %>%
rbind(c("3e: ANCOVA, baseline som kovariat", "-5.20 (-12.76, 2.37)", 0.14))
pred_df <- expand.grid(treat = c("P", "A"),
baseline = 10:30) %>%
mutate(pressure_pred = predict(ancova, newdata = .))
ggplot(filter(lang, minutter == 5), aes(x = baseline, y = pressure, colour = treat)) +
geom_point() +
geom_line(aes(y = pressure_pred), pred_df) +
labs(x = "Baseline pressure", y = "Pressure at 5 minuttes")
The plot doesn’t convey much information but serves only to show what we model: the treatment effect, which is the (constant!) vertical difference between the curves.
From now we look at all measurements and not just a specific point in time.
lang <- lang %>%
mutate(treatA_minutter = ifelse(treat == "A", minutter, 0),
minutter_fct = factor(minutter),
treatA_minutter_fct = factor(treatA_minutter))
xtabs(~ treat + treatA_minutter, data = filter(lang, minutter <= 5))
## treatA_minutter
## treat 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
## P 88 0 0 0 0 0 0 0 0 0 0
## A 7 7 7 7 7 7 7 6 7 6 5
model4 <- lme(pressure ~ minutter_fct + treatA_minutter_fct,
random = ~ 1 | patient,
na.action = na.exclude,
data = filter(lang, minutter <= 5))
summary(model4) # disregard the correlation matrix that just clutters everything
## Linear mixed-effects model fit by REML
## Data: filter(lang, minutter <= 5)
## AIC BIC logLik
## 699.3604 767.0182 -326.6802
##
## Random effects:
## Formula: ~1 | patient
## (Intercept) Residual
## StdDev: 7.12126 1.63761
##
## Fixed effects: pressure ~ minutter_fct + treatA_minutter_fct
## Value Std.Error DF t-value p-value
## (Intercept) 20.162500 1.8267819 125 11.037169 0.0000
## minutter_fct0.5 0.275554 0.7677241 125 0.358923 0.7203
## minutter_fct1 0.486665 0.7677241 125 0.633906 0.5273
## minutter_fct1.5 1.164443 0.7677241 125 1.516747 0.1319
## minutter_fct2 1.364443 0.7677241 125 1.777257 0.0780
## minutter_fct2.5 1.042221 0.7677241 125 1.357546 0.1771
## minutter_fct3 0.608887 0.7677241 125 0.793107 0.4292
## minutter_fct3.5 0.819998 0.7677241 125 1.068090 0.2875
## minutter_fct4 0.642221 0.7677241 125 0.836525 0.4045
## minutter_fct4.5 1.574489 1.0015045 125 1.572124 0.1184
## minutter_fct5 1.758823 1.1174624 125 1.573944 0.1180
## treatA_minutter_fct0.5 -0.244123 1.1523709 125 -0.211845 0.8326
## treatA_minutter_fct1 -1.012377 1.1523709 125 -0.878517 0.3813
## treatA_minutter_fct1.5 -1.604441 1.1523709 125 -1.392296 0.1663
## treatA_minutter_fct2 -2.761584 1.1523709 125 -2.396437 0.0180
## treatA_minutter_fct2.5 -3.353647 1.1523709 125 -2.910215 0.0043
## treatA_minutter_fct3 -3.906028 1.1523709 125 -3.389558 0.0009
## treatA_minutter_fct3.5 -4.721591 1.1824685 125 -3.992995 0.0001
## treatA_minutter_fct4 -5.339362 1.1523709 125 -4.633371 0.0000
## treatA_minutter_fct4.5 -7.324097 1.3462905 125 -5.440205 0.0000
## treatA_minutter_fct5 -8.423005 1.4685338 125 -5.735656 0.0000
## Correlation:
## (Intr) mn_0.5 mntt_1 mn_1.5 mntt_2 mn_2.5 mntt_3 mn_3.5
## minutter_fct0.5 -0.120
## minutter_fct1 -0.120 0.494
## minutter_fct1.5 -0.120 0.494 0.494
## minutter_fct2 -0.120 0.494 0.494 0.494
## minutter_fct2.5 -0.120 0.494 0.494 0.494 0.494
## minutter_fct3 -0.120 0.494 0.494 0.494 0.494 0.494
## minutter_fct3.5 -0.120 0.494 0.494 0.494 0.494 0.494 0.494
## minutter_fct4 -0.120 0.494 0.494 0.494 0.494 0.494 0.494 0.494
## minutter_fct4.5 -0.092 0.379 0.379 0.379 0.379 0.379 0.379 0.379
## minutter_fct5 -0.082 0.340 0.340 0.340 0.340 0.340 0.340 0.340
## treatA_minutter_fct0.5 0.000 -0.657 -0.320 -0.320 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct1 0.000 -0.320 -0.657 -0.320 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct1.5 0.000 -0.320 -0.320 -0.657 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct2 0.000 -0.320 -0.320 -0.320 -0.657 -0.320 -0.320 -0.320
## treatA_minutter_fct2.5 0.000 -0.320 -0.320 -0.320 -0.320 -0.657 -0.320 -0.320
## treatA_minutter_fct3 0.000 -0.320 -0.320 -0.320 -0.320 -0.320 -0.657 -0.320
## treatA_minutter_fct3.5 0.000 -0.312 -0.312 -0.312 -0.312 -0.312 -0.312 -0.640
## treatA_minutter_fct4 0.000 -0.320 -0.320 -0.320 -0.320 -0.320 -0.320 -0.320
## treatA_minutter_fct4.5 0.000 -0.274 -0.274 -0.274 -0.274 -0.274 -0.274 -0.274
## treatA_minutter_fct5 0.000 -0.251 -0.251 -0.251 -0.251 -0.251 -0.251 -0.251
## mntt_4 mn_4.5 mntt_5 tA__0. trA__1 tA__1. trA__2 tA__2.
## minutter_fct0.5
## minutter_fct1
## minutter_fct1.5
## minutter_fct2
## minutter_fct2.5
## minutter_fct3
## minutter_fct3.5
## minutter_fct4
## minutter_fct4.5 0.379
## minutter_fct5 0.340 0.297
## treatA_minutter_fct0.5 -0.320 -0.245 -0.220
## treatA_minutter_fct1 -0.320 -0.245 -0.220 0.487
## treatA_minutter_fct1.5 -0.320 -0.245 -0.220 0.487 0.487
## treatA_minutter_fct2 -0.320 -0.245 -0.220 0.487 0.487 0.487
## treatA_minutter_fct2.5 -0.320 -0.245 -0.220 0.487 0.487 0.487 0.487
## treatA_minutter_fct3 -0.320 -0.245 -0.220 0.487 0.487 0.487 0.487 0.487
## treatA_minutter_fct3.5 -0.312 -0.239 -0.214 0.475 0.475 0.475 0.475 0.475
## treatA_minutter_fct4 -0.657 -0.245 -0.220 0.487 0.487 0.487 0.487 0.487
## treatA_minutter_fct4.5 -0.274 -0.738 -0.215 0.417 0.417 0.417 0.417 0.417
## treatA_minutter_fct5 -0.251 -0.220 -0.756 0.382 0.382 0.382 0.382 0.382
## trA__3 tA__3. trA__4 tA__4.
## minutter_fct0.5
## minutter_fct1
## minutter_fct1.5
## minutter_fct2
## minutter_fct2.5
## minutter_fct3
## minutter_fct3.5
## minutter_fct4
## minutter_fct4.5
## minutter_fct5
## treatA_minutter_fct0.5
## treatA_minutter_fct1
## treatA_minutter_fct1.5
## treatA_minutter_fct2
## treatA_minutter_fct2.5
## treatA_minutter_fct3
## treatA_minutter_fct3.5 0.475
## treatA_minutter_fct4 0.487 0.475
## treatA_minutter_fct4.5 0.417 0.406 0.417
## treatA_minutter_fct5 0.382 0.371 0.382 0.352
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.96134903 -0.51474694 0.01765179 0.47638279 2.40697726
##
## Number of Observations: 161
## Number of Groups: 16
resids <- residuals(model4, level = 0) # population level
cov <- extract.lme.cov(model4)
L <- t(chol(cov)) # transposed cholesky defactorisation of the covariance matrix
plot_df <- tibble(scaled_residuals = solve(L) %*% resids, # solve(L) yields the inversed L
fitted_values = fitted(model4,level = 0))
ggplot(plot_df) +
geom_qq_line(aes(sample = scaled_residuals), colour = "grey60") +
geom_qq(aes(sample = scaled_residuals), size = 0.5)
ggplot(plot_df, aes(x = fitted_values, y = scaled_residuals)) +
geom_jitter(size = 0.5) +
geom_smooth()
qplot(x = fitted(model4,level=1), y = scale(resid(model4,level=1))) +
geom_hline(yintercept = 0)
qplot(sample = scale(resid(model4)), geom = "qq") +
geom_qq_line()
qplot(scale(resid(model4)))
Tungen lige i munden: minutter_fct* correspond to the deviation from the common mean of the placebo group, whereas the treatA_minutter_* corresponds to the difference between the placebo group and active arm. Thus, we can’t give a single estimate but must resort so timepoint-specific estimates because we didn’t make any assumptions about an overall trend from which the two groups could deviate. And we can see significant difference from the 2-minute mark onward.
intervals(model4)
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 16.5470719 20.1625000 23.7779281
## minutter_fct0.5 -1.2438673 0.2755540 1.7949753
## minutter_fct1 -1.0327561 0.4866651 2.0060864
## minutter_fct1.5 -0.3549784 1.1644429 2.6838642
## minutter_fct2 -0.1549784 1.3644429 2.8838642
## minutter_fct2.5 -0.4772006 1.0422207 2.5616420
## minutter_fct3 -0.9105339 0.6088874 2.1283086
## minutter_fct3.5 -0.6994228 0.8199985 2.3394197
## minutter_fct4 -0.8772006 0.6422207 2.1616420
## minutter_fct4.5 -0.4076123 1.5744894 3.5565911
## minutter_fct5 -0.4527737 1.7588231 3.9704198
## treatA_minutter_fct0.5 -2.5248085 -0.2441235 2.0365615
## treatA_minutter_fct1 -3.2930625 -1.0123774 1.2683076
## treatA_minutter_fct1.5 -3.8851260 -1.6044409 0.6762441
## treatA_minutter_fct2 -5.0422688 -2.7615838 -0.4808988
## treatA_minutter_fct2.5 -5.6343323 -3.3536473 -1.0729623
## treatA_minutter_fct3 -6.1867133 -3.9060282 -1.6253432
## treatA_minutter_fct3.5 -7.0618425 -4.7215906 -2.3813387
## treatA_minutter_fct4 -7.6200466 -5.3393616 -3.0586766
## treatA_minutter_fct4.5 -9.9885732 -7.3240971 -4.6596211
## treatA_minutter_fct5 -11.3294155 -8.4230049 -5.5165942
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: patient
## lower est. upper
## sd((Intercept)) 4.96722 7.12126 10.2094
##
## Within-group standard error:
## lower est. upper
## 1.446722 1.637610 1.853684
intervals(model4)$fixed[21, ] # estimates at specific time point
## lower est. upper
## -11.329415 -8.423005 -5.516594
The estimates are much different from the naive ones found earlier (3d) but more like those of 3e, however with much narrower confidence interval. This is because the current analysis included much more data, and this lends information to each specific time point.
We add the results to the table:
result <- result %>%
rbind(c("4b: CS-model med minutter som class", "-8.42 (-11.33, -5.53)", 0.0000))
There is a significant difference between the two groups overall.
model4_ml <- lme(pressure ~ minutter_fct + treatA_minutter_fct,
random = ~ 1 | patient,
na.action = na.exclude,
data = filter(lang, minutter <= 5),
method = "ML")
model4_red <- update(model4_ml, . ~ minutter_fct)
anova(model4_ml, model4_red)
## Model df AIC BIC logLik Test L.Ratio p-value
## model4_ml 1 23 723.1673 794.0396 -338.5837
## model4_red 2 13 775.6512 815.7095 -374.8256 1 vs 2 72.4839 <.0001
Models with lower AIC and BIC (but larger log-likelihood) are preferred, so modelling the difference between groups improves the model significantly. AIC and BIC account somewhat for the fact that a model with many parameters needs to be a lot better to be preferred.
plot_df <- model4$data %>%
mutate(pred4_fix = predict(model4,level=0), # fixed effects
pred4 = fitted(model4)) # patient-level effects
ggplot(plot_df, aes(x = minutter, y = pred4_fix, colour = treat)) +
geom_line() +
geom_point() +
coord_cartesian(ylim = c(0, NA))
ggplot(plot_df, aes(x = minutter, y = pred4, colour = patient)) +
geom_line() +
geom_point() +
facet_wrap(~ treat, nrow = 1)
Three things to notice:
(Intercept) of the StdDev parameter is close to 0).model5a <- lme(pressure ~ minutter_fct + treatA_minutter_fct,
random = ~1 | patient,
na.action = na.exclude,
correlation = corCAR1(0.9, ~ minutter | patient, FALSE),
data = filter(lang, minutter <= 5))
summary(model5a)
## Linear mixed-effects model fit by REML
## Data: filter(lang, minutter <= 5)
## AIC BIC logLik
## 611.0626 681.662 -281.5313
##
## Random effects:
## Formula: ~1 | patient
## (Intercept) Residual
## StdDev: 0.01361138 6.854021
##
## Correlation Structure: Continuous AR(1)
## Formula: ~minutter | patient
## Parameter estimate(s):
## Phi
## 0.9639969
## Fixed effects: pressure ~ minutter_fct + treatA_minutter_fct
## Value Std.Error DF t-value p-value
## (Intercept) 20.162500 1.7135087 125 11.766792 0.0000
## minutter_fct0.5 0.305343 0.4346212 125 0.702549 0.4836
## minutter_fct1 0.499881 0.6106497 125 0.818604 0.4146
## minutter_fct1.5 1.161386 0.7430574 125 1.562983 0.1206
## minutter_fct2 1.345410 0.8525002 125 1.578193 0.1170
## minutter_fct2.5 1.007502 0.9470442 125 1.063838 0.2895
## minutter_fct3 0.558767 1.0308598 125 0.542040 0.5888
## minutter_fct3.5 0.754757 1.1064440 125 0.682146 0.4964
## minutter_fct4 0.562132 1.1754375 125 0.478233 0.6333
## minutter_fct4.5 1.724982 1.3304142 125 1.296575 0.1972
## minutter_fct5 1.165922 1.5122688 125 0.770975 0.4422
## treatA_minutter_fct0.5 -0.312211 0.6553980 125 -0.476369 0.6346
## treatA_minutter_fct1 -1.042584 0.9184922 125 -1.135104 0.2585
## treatA_minutter_fct1.5 -1.597455 1.1148100 125 -1.432939 0.1544
## treatA_minutter_fct2 -2.718080 1.2757755 125 -2.130532 0.0351
## treatA_minutter_fct2.5 -3.274290 1.4137012 125 -2.316111 0.0222
## treatA_minutter_fct3 -3.791468 1.5349735 125 -2.470054 0.0149
## treatA_minutter_fct3.5 -4.334556 1.6496563 125 -2.627551 0.0097
## treatA_minutter_fct4 -5.156303 1.7415996 125 -2.960671 0.0037
## treatA_minutter_fct4.5 -7.795249 1.9049122 125 -4.092183 0.0001
## treatA_minutter_fct5 -7.802605 2.0976402 125 -3.719706 0.0003
## Correlation:
## (Intr) mn_0.5 mntt_1 mn_1.5 mntt_2 mn_2.5 mntt_3 mn_3.5
## minutter_fct0.5 -0.072
## minutter_fct1 -0.101 0.702
## minutter_fct1.5 -0.123 0.570 0.811
## minutter_fct2 -0.142 0.490 0.698 0.860
## minutter_fct2.5 -0.158 0.436 0.620 0.764 0.888
## minutter_fct3 -0.173 0.395 0.562 0.693 0.806 0.907
## minutter_fct3.5 -0.187 0.363 0.517 0.638 0.741 0.834 0.920
## minutter_fct4 -0.199 0.338 0.481 0.593 0.689 0.775 0.855 0.929
## minutter_fct4.5 -0.196 0.295 0.419 0.517 0.601 0.676 0.745 0.810
## minutter_fct5 -0.190 0.256 0.364 0.449 0.522 0.587 0.647 0.704
## treatA_minutter_fct0.5 0.000 -0.660 -0.461 -0.372 -0.318 -0.281 -0.254 -0.232
## treatA_minutter_fct1 0.000 -0.462 -0.658 -0.531 -0.454 -0.402 -0.362 -0.331
## treatA_minutter_fct1.5 0.000 -0.374 -0.532 -0.656 -0.562 -0.496 -0.448 -0.410
## treatA_minutter_fct2 0.000 -0.321 -0.457 -0.563 -0.655 -0.579 -0.522 -0.477
## treatA_minutter_fct2.5 0.000 -0.284 -0.405 -0.499 -0.580 -0.653 -0.589 -0.539
## treatA_minutter_fct3 0.000 -0.257 -0.366 -0.451 -0.525 -0.591 -0.651 -0.596
## treatA_minutter_fct3.5 0.000 -0.235 -0.334 -0.412 -0.479 -0.540 -0.595 -0.647
## treatA_minutter_fct4 0.000 -0.218 -0.311 -0.383 -0.446 -0.502 -0.553 -0.602
## treatA_minutter_fct4.5 0.000 -0.196 -0.279 -0.344 -0.400 -0.450 -0.497 -0.540
## treatA_minutter_fct5 0.000 -0.175 -0.249 -0.307 -0.357 -0.402 -0.443 -0.482
## mntt_4 mn_4.5 mntt_5 tA__0. trA__1 tA__1. trA__2 tA__2.
## minutter_fct0.5
## minutter_fct1
## minutter_fct1.5
## minutter_fct2
## minutter_fct2.5
## minutter_fct3
## minutter_fct3.5
## minutter_fct4
## minutter_fct4.5 0.872
## minutter_fct5 0.757 0.868
## treatA_minutter_fct0.5 -0.215 -0.186 -0.161
## treatA_minutter_fct1 -0.306 -0.266 -0.229 0.701
## treatA_minutter_fct1.5 -0.379 -0.328 -0.284 0.567 0.809
## treatA_minutter_fct2 -0.441 -0.383 -0.331 0.486 0.694 0.858
## treatA_minutter_fct2.5 -0.498 -0.432 -0.373 0.431 0.615 0.760 0.886
## treatA_minutter_fct3 -0.551 -0.478 -0.413 0.390 0.556 0.687 0.801 0.904
## treatA_minutter_fct3.5 -0.598 -0.519 -0.448 0.356 0.508 0.628 0.732 0.826
## treatA_minutter_fct4 -0.648 -0.562 -0.486 0.331 0.472 0.584 0.681 0.768
## treatA_minutter_fct4.5 -0.582 -0.672 -0.580 0.297 0.424 0.524 0.611 0.690
## treatA_minutter_fct5 -0.519 -0.599 -0.695 0.265 0.378 0.467 0.545 0.615
## trA__3 tA__3. trA__4 tA__4.
## minutter_fct0.5
## minutter_fct1
## minutter_fct1.5
## minutter_fct2
## minutter_fct2.5
## minutter_fct3
## minutter_fct3.5
## minutter_fct4
## minutter_fct4.5
## minutter_fct5
## treatA_minutter_fct0.5
## treatA_minutter_fct1
## treatA_minutter_fct1.5
## treatA_minutter_fct2
## treatA_minutter_fct2.5
## treatA_minutter_fct3
## treatA_minutter_fct3.5 0.914
## treatA_minutter_fct4 0.850 0.923
## treatA_minutter_fct4.5 0.763 0.829 0.898
## treatA_minutter_fct5 0.680 0.739 0.800 0.892
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5330963 -0.9428729 -0.3393333 0.7592998 2.0852057
##
## Number of Observations: 161
## Number of Groups: 16
intervals(model5a, which = "fixed")$fixed[21, ] # note we must specify which
## lower est. upper
## -11.954095 -7.802605 -3.651115
model5a_ml <- update(model5a, method = "ML")
model5a_red <- update(model5a, . ~ minutter_fct, method = "ML")
anova(model5a_ml, model5a_red) # still significantly better model *with* group differences
## Model df AIC BIC logLik Test L.Ratio p-value
## model5a_ml 1 24 620.8233 694.7770 -286.4117
## model5a_red 2 14 622.6560 665.7957 -297.3280 1 vs 2 21.83272 0.016
result <- result %>%
rbind(c("5: SP(POW)-model med minutter som class", "-7.80 (-12.00, -3.61)", 0.0003))
compute_residuals <- function(model) {
large_resids <- residuals(model, level = 0)
cov <- extract.lme.cov(model, lang)
L <- t(chol(cov))
model$data %>%
mutate(large_resids_scaled = as.numeric(solve(L) %*% large_resids),
conditional_resids_scaled = as.numeric(scale(resid(model))),
fitted_values = fitted(model, level = 0))
}
ggplot(compute_residuals(model5a), aes(sample = large_resids_scaled)) +
geom_qq() +
geom_qq_line()
ggplot(compute_residuals(model5a), aes(x = fitted_values, y= large_resids_scaled)) +
geom_point() +
geom_smooth()
ggplot(compute_residuals(model5a), aes(sample = conditional_resids_scaled)) +
geom_qq() +
geom_qq_line()
plot_df <- model5a$data %>%
mutate(pred5_fix = predict(model5a, level = 0), # fixed effects
pred5 = fitted(model5a)) # patient-level effects
ggplot(plot_df, aes(x = minutter, y = pred5_fix, colour = treat)) +
geom_line() +
geom_point()
# This plot shows that the predictions for all patients are almost the same unlike model4
ggplot(plot_df, aes(x = minutter, y = pred5, colour = patient)) +
geom_line() +
geom_point() +
facet_wrap(~ treat, nrow = 1)
model6a <- lme(pressure ~ minutter + treat:minutter,
random = ~ 1 | patient,
data = lang,
correlation = corCAR1(0.9, ~ minutter | patient, FALSE),
na.action = na.omit)
summary(model6a)
## Linear mixed-effects model fit by REML
## Data: lang
## AIC BIC logLik
## 601.4561 619.8316 -294.728
##
## Random effects:
## Formula: ~1 | patient
## (Intercept) Residual
## StdDev: 4.626822 5.02678
##
## Correlation Structure: Continuous AR(1)
## Formula: ~minutter | patient
## Parameter estimate(s):
## Phi
## 0.9342032
## Fixed effects: pressure ~ minutter + treat:minutter
## Value Std.Error DF t-value p-value
## (Intercept) 20.234116 1.7075505 143 11.849790 0.0000
## minutter 0.204885 0.2698103 143 0.759366 0.4489
## minutter:treatA -1.515154 0.3882902 143 -3.902118 0.0001
## Correlation:
## (Intr) minttr
## minutter -0.201
## minutter:treatA 0.002 -0.667
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.32425075 -0.63786627 -0.04937682 0.64238165 1.97106796
##
## Number of Observations: 161
## Number of Groups: 16
ggplot(compute_residuals(model6a), aes(sample = large_resids_scaled)) +
geom_qq() +
geom_qq_line()
ggplot(compute_residuals(model6a), aes(x = fitted_values, y= large_resids_scaled)) +
geom_point() +
geom_smooth()
ggplot(compute_residuals(model6a), aes(sample = conditional_resids_scaled)) +
geom_qq() +
geom_qq_line()
plot_df <- model6a$data %>%
mutate(pred_fix = predict(model6a, level = 0), # fixed effects
pred = fitted(model6a)) # patient-level effects
ggplot(plot_df, aes(x = minutter, y = pred_fix, colour = treat)) +
geom_line() +
geom_point()
ggplot(plot_df, aes(x = minutter, y = pred, colour = patient)) +
geom_line() +
geom_point() +
facet_wrap(~ treat, nrow = 1)
Compare with model 5A
model6a_ml = update(model6a, method = "ML")
anova(model5a_ml, model6a_ml)
## Model df AIC BIC logLik Test L.Ratio p-value
## model5a_ml 1 24 620.8233 694.7770 -286.4117
## model6a_ml 2 6 602.6345 621.1229 -295.3173 1 vs 2 17.8112 0.4682
# The estimated effect at 5 minutes (because linear we multiply by 5)
intervals(model6a, which="fixed")$fixed[3,] * 5
## lower est. upper
## -11.413423 -7.575772 -3.738121
result <- result %>%
rbind(c("6a: SP(POW)-model med lineær effekt af tid", "-7.58 (-11.41, -3.74)", 0.0001))
model6b <- lme(pressure ~ minutter + treat:minutter,
random = ~ 1 + minutter | patient,
na.action = na.omit,
data = lang)
summary(model6b)
## Linear mixed-effects model fit by REML
## Data: lang
## AIC BIC logLik
## 631.9417 653.3798 -308.9708
##
## Random effects:
## Formula: ~1 + minutter | patient
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 7.4684951 (Intr)
## minutter 0.7919332 -0.289
## Residual 1.0937369
##
## Fixed effects: pressure ~ minutter + treat:minutter
## Value Std.Error DF t-value p-value
## (Intercept) 20.761744 1.8739862 143 11.078920 0.0000
## minutter 0.172533 0.2719200 143 0.634499 0.5268
## minutter:treatA -1.564592 0.3981502 143 -3.929652 0.0001
## Correlation:
## (Intr) minttr
## minutter -0.226
## minutter:treatA 0.002 -0.648
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.02054204 -0.48388203 -0.03573534 0.49987358 3.00947584
##
## Number of Observations: 161
## Number of Groups: 16
ggplot(compute_residuals(model6b), aes(sample = large_resids_scaled)) +
geom_qq() +
geom_qq_line()
ggplot(compute_residuals(model6b), aes(x = fitted_values, y= large_resids_scaled)) +
geom_point() +
geom_smooth()
ggplot(compute_residuals(model6b), aes(sample = conditional_resids_scaled)) +
geom_qq() +
geom_qq_line()
ggplot(compute_residuals(model6b), aes(x = minutter, y = conditional_resids_scaled, colour = patient)) +
geom_line() +
geom_point() +
facet_wrap(~ treat, nrow = 1)
plot_df <- model6b$data %>%
mutate(pred_fix = predict(model6b, level = 0), # fixed effects
pred = fitted(model6b)) # patient-level effects
ggplot(plot_df, aes(x = minutter, y = pred_fix, colour = treat)) +
geom_line() +
geom_point()
ggplot(plot_df, aes(x = minutter, y = pred, colour = patient)) +
geom_line() +
geom_point() +
facet_wrap(~ treat, nrow = 1)
# Again, linear model so we just multiply the coefficient by 5
intervals(model6b, which = "fixed")$fixed[3, ] * 5
## lower est. upper
## -11.758061 -7.822959 -3.887857
result <- result %>%
rbind(c("6b: Random-effects regression", "-7.82 (-11.76, -3.89)", 0.0001))
Models 6a and 6b have the same mean-value structure, so we can compare them directly. The AR(1) model has a lower AIC so we prefer that as it fits the data better.
anova(model6a, model6b)
## Model df AIC BIC logLik Test L.Ratio p-value
## model6a 1 6 601.4561 619.8316 -294.7280
## model6b 2 7 631.9417 653.3798 -308.9708 1 vs 2 28.48562 <.0001
Including the outlier doesn’t change the overall conclusion (active-arm patients have a steeper trend) and the trend is even amplified. This makes sense if we consider the data of that outlier.
model7 = lme(pressure ~ minutter + treat:minutter,
random = ~ 1 | patient, data = lang5,
correlation = corCAR1(0.9, ~ minutter | patient, FALSE),
na.action = na.omit)
summary(model7)
## Linear mixed-effects model fit by REML
## Data: lang5
## AIC BIC logLik
## 661.4622 680.2059 -324.7311
##
## Random effects:
## Formula: ~1 | patient
## (Intercept) Residual
## StdDev: 5.057275 4.41461
##
## Correlation Structure: Continuous AR(1)
## Formula: ~minutter | patient
## Parameter estimate(s):
## Phi
## 0.9011342
## Fixed effects: pressure ~ minutter + treat:minutter
## Value Std.Error DF t-value p-value
## (Intercept) 20.362616 1.6271077 152 12.514608 0.0000
## minutter 0.192781 0.2812642 152 0.685409 0.4941
## minutter:treatA -1.614489 0.3907038 152 -4.132258 0.0001
## Correlation:
## (Intr) minttr
## minutter -0.209
## minutter:treatA 0.002 -0.689
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.42325828 -0.57640556 0.02810233 0.55744791 1.91974344
##
## Number of Observations: 171
## Number of Groups: 17
intervals(model7)$fixed[3, ] * 5
## lower est. upper
## -11.932000 -8.072444 -4.212889
result <- result %>%
rbind(c("7: som 6a, med person nr. 13", "-8.07 (-11.93, -4.21)", 0.0001))
result
## # A tibble: 8 x 3
## metode behandlingseffekt_5_minutt… p_value
## <chr> <chr> <chr>
## 1 3d: T-test efter 5 minutter -0.27 (-7.57, 7.01) 0.92
## 2 3e: T-test på ændring 0-5 minutter -7.61 (-12.88, -2.33) 0.015
## 3 3e: ANCOVA, baseline som kovariat -5.20 (-12.76, 2.37) 0.14
## 4 4b: CS-model med minutter som class -8.42 (-11.33, -5.53) 0
## 5 5: SP(POW)-model med minutter som class -7.80 (-12.00, -3.61) 3e-04
## 6 6a: SP(POW)-model med lineær effekt af tid -7.58 (-11.41, -3.74) 1e-04
## 7 6b: Random-effects regression -7.82 (-11.76, -3.89) 1e-04
## 8 7: som 6a, med person nr. 13 -8.07 (-11.93, -4.21) 1e-04
# diagnostics as for model 6a